exp: Padé denominator solve via LinearSolve default algorithm choice#239
Merged
ChrisRackauckas merged 5 commits intoJul 5, 2026
Merged
Conversation
Contributor
Author
|
Benchmarked honestly: as written, this is a performance regression in the hot regime — do not merge yet.
Slower up to ~1.7× for n ≤ 60 (the Krylov- Next revision here: thread a cached |
ExpMethodHigham2005Base solved (V - U) X = (V + U) with raw LAPACK.gesv! calls; route both sites (small-norm and squaring branches) through LinearSolve's default algorithm choice with the batched matrix RHS instead. This picks the best factorization per size/type (e.g. RecursiveFactorization for small matrices when loaded), and fails with a retcode rather than a LAPACKException on pathological input. Verified exp ≈ Base.exp for n = 5, 40, 120 at small and large norms (both code paths) against LinearSolve 4.0.0. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh
The per-call solve(LinearProblem(temp, X)) added init, matrix copies, and a fresh factorization allocation on every exponential! call, which measured slower than the raw LAPACK.gesv! it replaced for n <= 60. Build a LinearSolve.init workspace once in alloc_mem with workspace-owned n-by-n buffers and alias = LinearAliasSpecifier( alias_A = true, alias_b = true), so each call just refills the buffers and solve! refactorizes in place (lu!) with no O(n^2) copy (LinearSolve >= 4.1). Both Pade branches (small-norm and scaling-and-squaring) route through the shared helper, which restores gesv!'s SingularException contract on a failed factorization. Best-of-N per-call times (Julia 1.12, OpenBLAS, this machine), old gesv! vs cached LinearSolve: ~1.3-1.8x slower for n <= 30 (fixed solve!/property overhead of a few microseconds), parity at n = 60, and 1.2-1.5x faster for n >= 120. Allocations on cache reuse stay O(n) (e.g. 1.4 KB at n = 60). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Drop the pinned LUFactorization in alloc_mem and let LinearSolve.init size-select the algorithm (GenericLUFactorization for n <= 10, tuned LAPACK/RecursiveFactorization choices above). The tiny-matrix generic LU refactorizes in place with a cached pivot vector, making the warm solve! allocation-free and at parity with (n = 5) or faster than (n = 10) the raw LAPACK.gesv! this path originally used. Best-of-N end-to-end exponential! times against the old gesv! code (Julia 1.12, OpenBLAS single-threaded, this machine) are within ~3% at all of n = 5, 10, 30, 60, 120, 250 for scales 0.3 and 4.0, with fewer per-call allocation bytes at n <= 10. One behavioral note: an exactly singular Pade denominator is now resolved by the default algorithm's QR safety fallback instead of throwing, matching generic \ semantics. The denominator is nonsingular by construction for the branch norm bounds, so the retcode guard only fires if the fallback itself fails. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ChrisRackauckas-Claude
force-pushed
the
exp-linearsolve
branch
from
a312efa to
70b6b93
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stacked on #236 (branches from it; shares the LinearSolve v4 dependency — merge #236 first).
First step of moving the package's remaining linear solves onto LinearSolve:
ExpMethodHigham2005Base's Padé denominator solve(V−U) X = (V+U)(both the small-norm and squaring branches) now goes through LinearSolve's default algorithm choice with the v4 batched matrix RHS, replacing the rawLAPACK.gesv!calls. Benefits: per-size/type algorithm selection (e.g. RecursiveFactorization when loaded), retcode-based failure instead ofLAPACKException, consistency withphi!.Verified locally against LinearSolve 4.0.0:
exponential! ≈ Base.expfor n = 5, 40, 120 at small and large norms (exercising both converted sites).Scoped small per review preference. Follow-ups, in order: (1)
ExpMethodHigham2005(exp_noalloc.jl) with a workspace-embedded LinearSolve cache once the upstream in-place dense-LU option lands (SciML/LinearSolve.jl#1074 discussion); (2) benchmark whether a cachedinit/solve!here is worth the extraalloc_memplumbing. Deliberately not convertingexp_generic.jl(exotic-eltype path; generic LU is already optimal there) or the SMatrix\(container-preserving).CI note: needs LinearSolve 4.0.0 registered in General (in flight) to resolve.
🤖 Generated with Claude Code